To install sf from CRAN:
install.packages("sf")
or the development version from GitHub:
devtools::install_github("edzer/sfr")
then load:
library(sf)
NB MacOSX and LINUX users need to install a number of geospatial libraries (GEOS, GDAL, and proj.4).
The tidyverse package also needs to be loaded:
library(tidyverse)
A GeoJSON of Greater Manchester’s wards was created from a vector boundary file available from ONS’s Open Geography Portal. The GeoJSON is projected in British National Grid (EPSG:27700) and originally derives from the Ordnance Survey.
Point data are supplied by data.police.uk and represent incidents of anti-social behaviour and crime recorded by the Greater Manchester Police during February 2017. The incidents are supplied with latitude and longitude coordinates which have undergone an anonymisation process.
Reading spatial data (polygons)
bdy <- st_read("data/wards.geojson")
Reading layer `OGRGeoJSON' from data source `/Users/hcpartridge/Desktop/rumgroup/data/wards.geojson' using driver `GeoJSON'
converted into: POLYGON
Simple feature collection with 215 features and 12 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 351664 ymin: 381168.6 xmax: 406087.5 ymax: 421039.8
epsg (SRID): 27700
proj4string: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
bdy
Simple feature collection with 215 features and 12 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 351664 ymin: 381168.6 xmax: 406087.5 ymax: 421039.8
epsg (SRID): 27700
proj4string: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
objectid wd16cd wd16nm wd16nmw lad16cd lad16nm bng_e bng_n long
1 772 E05000851 Ince <NA> E08000010 Wigan 359959 405278 -2.60570
2 773 E05000852 Leigh East <NA> E08000010 Wigan 367298 400695 -2.49447
3 774 E05000853 Leigh South <NA> E08000010 Wigan 366260 398660 -2.50990
4 775 E05000854 Leigh West <NA> E08000010 Wigan 363428 400762 -2.55282
5 776 E05000855 Lowton East <NA> E08000010 Wigan 362642 397699 -2.56431
6 777 E05000856 Orrell <NA> E08000010 Wigan 353321 404083 -2.70568
7 778 E05000857 Pemberton <NA> E08000010 Wigan 355101 405469 -2.67903
8 779 E05000858 Shevington with Lower Ground <NA> E08000010 Wigan 354692 408504 -2.68564
9 780 E05000859 Standish with Langtree <NA> E08000010 Wigan 355996 410264 -2.66620
10 781 E05000860 Tyldesley <NA> E08000010 Wigan 369533 402337 -2.46094
11 782 E05000861 Wigan Central <NA> E08000010 Wigan 358550 406484 -2.62712
12 783 E05000862 Wigan West <NA> E08000010 Wigan 356532 406901 -2.65763
13 784 E05000863 Winstanley <NA> E08000010 Wigan 355407 402755 -2.67402
14 785 E05000864 Worsley Mesnes <NA> E08000010 Wigan 357808 403626 -2.63793
15 571 E05000650 Astley Bridge <NA> E08000001 Bolton 370670 412906 -2.44479
16 572 E05000651 Bradshaw <NA> E08000001 Bolton 374590 413058 -2.38555
17 573 E05000652 Breightmet <NA> E08000001 Bolton 374437 409641 -2.38758
18 574 E05000653 Bromley Cross <NA> E08000001 Bolton 372113 414150 -2.42309
19 575 E05000654 Crompton <NA> E08000001 Bolton 371741 410493 -2.42838
20 576 E05000655 Farnworth <NA> E08000001 Bolton 373564 406162 -2.40047
21 577 E05000656 Great Lever <NA> E08000001 Bolton 371959 408204 -2.42488
22 578 E05000657 Halliwell <NA> E08000001 Bolton 370828 409676 -2.44210
23 579 E05000658 Harper Green <NA> E08000001 Bolton 371512 406080 -2.43144
24 580 E05000659 Heaton and Lostock <NA> E08000001 Bolton 367961 409240 -2.48536
25 581 E05000660 Horwich and Blackrod <NA> E08000001 Bolton 362140 410545 -2.57343
26 582 E05000661 Horwich North East <NA> E08000001 Bolton 365750 412038 -2.51906
27 583 E05000662 Hulton <NA> E08000001 Bolton 369418 405943 -2.46303
28 584 E05000663 Kearsley <NA> E08000001 Bolton 375564 404886 -2.37019
29 585 E05000664 Little Lever and Darcy Lever <NA> E08000001 Bolton 375274 407414 -2.37476
30 586 E05000665 Rumworth <NA> E08000001 Bolton 370267 408025 -2.45041
31 587 E05000666 Smithills <NA> E08000001 Bolton 369230 411951 -2.46646
32 588 E05000667 Tonge with the Haulgh <NA> E08000001 Bolton 372953 409770 -2.41001
33 589 E05000668 Westhoughton North and Chew Moor <NA> E08000001 Bolton 365620 406875 -2.52045
34 590 E05000669 Westhoughton South <NA> E08000001 Bolton 364395 404864 -2.53871
35 591 E05000670 Besses <NA> E08000002 Bury 381771 405926 -2.27659
36 592 E05000671 Church <NA> E08000002 Bury 378564 410535 -2.32531
37 593 E05000672 East <NA> E08000002 Bury 382701 411330 -2.26286
38 594 E05000673 Elton <NA> E08000002 Bury 379386 412322 -2.31301
39 595 E05000674 Holyrood <NA> E08000002 Bury 383274 405503 -2.25388
40 596 E05000675 Moorside <NA> E08000002 Bury 381103 412374 -2.28706
41 597 E05000676 North Manor <NA> E08000002 Bury 378620 414508 -2.32474
42 598 E05000677 Pilkington Park <NA> E08000002 Bury 379384 405697 -2.31260
43 599 E05000678 Radcliffe East <NA> E08000002 Bury 378655 408195 -2.32378
44 600 E05000679 Radcliffe North <NA> E08000002 Bury 376513 409386 -2.35621
45 601 E05000680 Radcliffe West <NA> E08000002 Bury 377876 406692 -2.33543
46 602 E05000681 Ramsbottom <NA> E08000002 Bury 378147 417005 -2.33207
47 603 E05000682 Redvales <NA> E08000002 Bury 380506 409255 -2.29589
48 604 E05000683 St Mary's <NA> E08000002 Bury 380107 403393 -2.30154
49 605 E05000684 Sedgley <NA> E08000002 Bury 382536 403044 -2.26488
50 606 E05000685 Tottington <NA> E08000002 Bury 377233 412927 -2.34559
51 607 E05000686 Unsworth <NA> E08000002 Bury 382436 407933 -2.26667
52 608 E05000687 Ancoats and Clayton <NA> E08000003 Manchester 387636 398888 -2.18779
53 609 E05000688 Ardwick <NA> E08000003 Manchester 385839 396661 -2.21477
54 610 E05000689 Baguley <NA> E08000003 Manchester 381171 388024 -2.28457
55 611 E05000690 Bradford <NA> E08000003 Manchester 387172 398057 -2.19475
56 612 E05000691 Brooklands <NA> E08000003 Manchester 380793 389733 -2.29036
57 613 E05000692 Burnage <NA> E08000003 Manchester 386408 392483 -2.20602
58 614 E05000693 Charlestown <NA> E08000003 Manchester 387385 403256 -2.19175
59 615 E05000694 Cheetham <NA> E08000003 Manchester 384408 400457 -2.23652
60 616 E05000695 Chorlton <NA> E08000003 Manchester 381131 393791 -2.28552
61 617 E05000696 Chorlton Park <NA> E08000003 Manchester 382585 392536 -2.26356
62 618 E05000697 City Centre <NA> E08000003 Manchester 383942 398119 -2.24342
63 619 E05000698 Crumpsall <NA> E08000003 Manchester 384444 402511 -2.23608
64 620 E05000699 Didsbury East <NA> E08000003 Manchester 384893 390817 -2.22874
65 621 E05000700 Didsbury West <NA> E08000003 Manchester 383755 391191 -2.24588
66 622 E05000701 Fallowfield <NA> E08000003 Manchester 384495 394300 -2.23490
67 623 E05000702 Gorton North <NA> E08000003 Manchester 388555 396421 -2.17385
68 624 E05000703 Gorton South <NA> E08000003 Manchester 388320 395173 -2.17734
69 625 E05000704 Harpurhey <NA> E08000003 Manchester 386364 401341 -2.20707
70 626 E05000705 Higher Blackley <NA> E08000003 Manchester 384361 404086 -2.23741
71 627 E05000706 Hulme <NA> E08000003 Manchester 383630 396698 -2.24805
72 628 E05000707 Levenshulme <NA> E08000003 Manchester 387196 393796 -2.19421
73 629 E05000708 Longsight <NA> E08000003 Manchester 387026 395294 -2.19684
74 630 E05000709 Miles Platting and Newton Heath <NA> E08000003 Manchester 387660 400026 -2.18748
75 631 E05000710 Moss Side <NA> E08000003 Manchester 383884 395430 -2.24416
76 632 E05000711 Moston <NA> E08000003 Manchester 387958 401864 -2.18305
lat st_areasha st_lengths geometry
1 53.54262 4780173 11442.625 POLYGON((359198.460435328 4...
2 53.50194 4972092 11056.203 POLYGON((367700.613520717 4...
3 53.48358 6521214 12377.967 POLYGON((364489.980933077 3...
4 53.50228 7315344 16159.924 POLYGON((363107.910792193 3...
5 53.47470 11185415 17990.730 POLYGON((361792.312268886 3...
6 53.53133 8884317 19650.851 POLYGON((353418.76434581 40...
7 53.54394 6239016 12622.229 POLYGON((355533.978524123 4...
8 53.57118 8963755 17516.213 POLYGON((353418.76434581 40...
9 53.58711 10401896 16021.993 POLYGON((356796.803981366 4...
10 53.51683 5559650 12889.007 POLYGON((368822.197554502 4...
11 53.55335 3525250 9567.055 POLYGON((357695.08368452 40...
12 53.55693 2708373 9178.338 POLYGON((356255.095937935 4...
13 53.51957 6183791 10253.866 POLYGON((355533.978524123 4...
14 53.52760 5930122 11368.272 POLYGON((358614.66424552 40...
15 53.61189 6601953 14566.948 POLYGON((370549.775116555 4...
16 53.61346 9283676 16445.880 POLYGON((375026.364570552 4...
17 53.58274 3704010 9695.609 POLYGON((373586.218118696 4...
18 53.62315 6975068 15981.476 POLYGON((373169.622402638 4...
19 53.59026 3483490 12335.029 POLYGON((372158.45537916 40...
20 53.55143 4355269 11771.147 POLYGON((374784.205484859 4...
21 53.56970 4036982 11531.263 POLYGON((370802.588433939 4...
22 53.58287 2820871 7928.388 POLYGON((369793.9063342 409...
23 53.55058 3158073 9727.627 POLYGON((370505.503263172 4...
24 53.57878 8065909 17328.300 POLYGON((366456.699782229 4...
25 53.59012 14745198 24960.144 POLYGON((363194.780215914 4...
26 53.60379 10174487 18821.484 POLYGON((366456.699782229 4...
27 53.54923 10475447 14821.862 POLYGON((368584.593338779 4...
28 53.54006 7692125 14032.845 POLYGON((376604.682156536 4...
29 53.56276 5163985 13467.401 POLYGON((374784.205484859 4...
30 53.56800 1937048 6183.851 POLYGON((370505.503263172 4...
31 53.60322 10977728 20306.973 POLYGON((366341.544716241 4...
32 53.58383 3078336 8853.175 POLYGON((372158.45537916 40...
33 53.55738 16751948 24438.248 POLYGON((363114.143528574 4...
34 53.53922 6313599 13102.652 POLYGON((363400.546355882 4...
35 53.54966 2544591 7588.175 POLYGON((380940.289808304 4...
36 53.59096 4036063 12604.877 POLYGON((375857.068369663 4...
37 53.59826 4195185 13450.955 POLYGON((379584.387973824 4...
38 53.60706 3316545 9011.894 POLYGON((378920.194781534 4...
39 53.54591 4460157 12444.692 POLYGON((382464.27927018 40...
40 53.60759 4703167 10772.104 POLYGON((381389.068949916 4...
41 53.62667 9789443 23196.083 POLYGON((375026.364570552 4...
42 53.54751 7217589 15340.667 POLYGON((380774.396924753 4...
43 53.56993 6470312 14319.963 POLYGON((379233.131869085 4...
44 53.58055 7053574 13177.442 POLYGON((375332.182010675 4...
45 53.55639 4896161 12607.624 POLYGON((377289.29754458 40...
46 53.64910 13489427 21951.402 POLYGON((375569.530811583 4...
47 53.57953 3676876 11852.667 POLYGON((381691.379115022 4...
48 53.52683 5014610 11047.551 POLYGON((380940.289808304 4...
49 53.52378 2110743 8605.973 POLYGON((381671.176295511 4...
50 53.61241 7792259 17558.920 POLYGON((375857.068369663 4...
51 53.56772 8696448 16181.589 POLYGON((381465.482109123 4...
52 53.48657 4434567 15297.173 POLYGON((388979.59511252 39...
53 53.46651 4090129 8290.054 POLYGON((385048.700403615 3...
54 53.38873 3819630 9646.251 POLYGON((382222.444481428 3...
55 53.47909 5160234 12321.011 POLYGON((389082.893464996 3...
56 53.40408 4166786 13208.062 POLYGON((379868.483208293 3...
57 53.42897 2324455 7870.979 POLYGON((385852.697738632 3...
58 53.52583 3572634 8915.117 POLYGON((388814.218263555 4...
59 53.50059 4703758 10305.930 POLYGON((385383.442175553 3...
60 53.44056 2495801 8555.880 POLYGON((382069.783909148 3...
61 53.42933 4798682 11399.285 POLYGON((382951.464352811 3...
62 53.47956 3102193 8569.712 POLYGON((384111.586234225 3...
63 53.51905 3186206 10323.461 POLYGON((383059.682257273 4...
64 53.41395 3542769 10341.897 POLYGON((384160.163612966 3...
65 53.41728 2823282 10970.634 POLYGON((384132.27552782 39...
66 53.44525 2265053 8571.091 POLYGON((382951.464352811 3...
67 53.46442 3520488 9673.390 POLYGON((390196.360627751 3...
68 53.45319 3315488 11618.528 POLYGON((387559.543108503 3...
69 53.50859 3815368 9202.353 POLYGON((385383.442175553 3...
70 53.53320 7040635 13779.614 POLYGON((382464.27927018 40...
71 53.46678 2396423 7806.111 POLYGON((383254.103669637 3...
72 53.44079 2356776 9676.619 POLYGON((386989.75973683 39...
73 53.45425 1498778 6469.738 POLYGON((385590.186839664 3...
74 53.49680 4353329 12542.810 POLYGON((384894.453135228 3...
75 53.45539 1825331 7563.015 POLYGON((383254.103669637 3...
76 53.51333 3093419 10756.913 POLYGON((388794.32215659 40...
[ reached getOption("max.print") -- omitted 139 rows ]
Reading data with coordinates (points)
pts <- read_csv("data/2017-03-greater-manchester-street.csv")
pts <- st_as_sf(pts, coords = c("Longitude", "Latitude"))
pts
Simple feature collection with 35094 features and 10 fields
geometry type: POINT
dimension: XY
bbox: xmin: -2.716562 ymin: 53.34144 xmax: -1.96473 ymax: 53.66877
epsg (SRID): NA
proj4string: NA
Writing spatial data
st_write(bdy, "boundaries.shp")
Attribute table (dataframe) AND geometry (list-column) with coordinates, CRS, bbox
st_geometry(bdy) # print geometry
Geometry set for 215 features
geometry type: POLYGON
dimension: XY
bbox: xmin: 351664 ymin: 381168.6 xmax: 406087.5 ymax: 421039.8
epsg (SRID): 27700
proj4string: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
First 5 geometries:
Rename variables
glimpse(bdy)
Observations: 215
Variables: 4
$ ward <fctr> Ince, Leigh East, Leigh South, Leigh West, Lowton East, Orrell, Pemberton, Shevington wi...
$ census_code <fctr> E05000851, E05000852, E05000853, E05000854, E05000855, E05000856, E05000857, E05000858, ...
$ borough <fctr> Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wiga...
$ geometry <simple_feature> POLYGON((359198.460435328 4..., POLYGON((367700.613520717 4..., POLYGON((36448...
Count frequency of wards by borough
bdy %>%
group_by(borough) %>%
count() %>%
arrange(desc(n)) %>% # sort in descending order
st_set_geometry(., NULL) # hide geometry
Select features
chorlton <- bdy %>%
filter(ward == "Chorlton")
plot(st_geometry(bdy))
plot(st_geometry(chorlton), col = "red", add = TRUE)
Using sf functions to add geometry column to dplyr chain
bdy <- bdy %>%
mutate(area = st_area(.)) # returns the area of a feature
bdy %>%
select(ward, area) %>%
arrange(desc(area)) %>%
slice(1:10) %>%
st_set_geometry(., NULL)
Check and assign CRS
st_crs(pts)
$epsg
[1] NA
$proj4string
[1] NA
attr(,"class")
[1] "crs"
pts <- st_set_crs(pts, 4326) # assign Lat/Long (epsg:4326)
st_crs(pts)
$epsg
[1] 4326
$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"
attr(,"class")
[1] "crs"
Reproject CRS
bdy_WGS84 <- st_transform(bdy, 4326)
st_crs(bdy_WGS84)
$epsg
[1] 4326
$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"
attr(,"class")
[1] "crs"
Convert to and from sp objects
bdy_sp <- as(bdy, 'Spatial')
class(bdy_sp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
bdy_sf <- st_as_sf(bdy_sp)
class(bdy_sf)
Buffer features
buffer <- chorlton %>%
st_buffer(dist = 1000)
plot(st_geometry(buffer))
plot(st_geometry(chorlton), col = "red", add = TRUE)
Buffer and intersect
pts_sub <- bdy %>%
filter(ward == "Chorlton") %>%
st_buffer(dist = 1000) %>%
st_intersection(st_transform(pts, 27700)) # reproject pts to BNG (epsg:27700)
plot(st_geometry(buffer))
plot(st_geometry(chorlton), col = "red", add = TRUE)
plot(st_geometry(pts_sub), col = "black", add = TRUE)
pts_sub %>%
group_by(`Crime type`) %>%
count() %>%
arrange(desc(n)) %>%
st_set_geometry(., NULL)
Points in polygon
pts %>%
filter(`Crime type` == "Vehicle crime") %>%
st_join(bdy_WGS84, ., left = FALSE) %>%
count(ward) %>%
arrange(desc(n)) %>%
st_set_geometry(., NULL)
bdy_pts <- pts %>%
filter(`Crime type` == "Vehicle crime") %>%
st_join(bdy_WGS84, ., left = FALSE) %>%
count(ward)
Base plots
plot(bdy_pts) # plots small multiples if dataframe has several attributes
plot(bdy_pts["n"]) # select the appropriate attribute to plot a single map
library(RColorBrewer) ; library(classInt)
pal <- brewer.pal(5, "RdPu")
classes <- classIntervals(bdy_pts$n, n=5, style="pretty")$brks
plot(bdy_pts["n"],
col = pal[findInterval(bdy_pts$n, classes, all.inside=TRUE)],
main = "Vehicle crime in Greater Manchester\nMarch 2017", axes = F)
legend("bottomright", legend = paste("<", round(classes[-1])), fill = pal, cex = 0.7)
devtools::install_github("tidyverse/ggplot2") # NB need development version for geom_sf()
ggplot(bdy_pts) +
geom_sf(aes(fill = n)) +
scale_fill_gradientn('Frequency', colours=RColorBrewer::brewer.pal(5,"RdPu"),
breaks = scales::pretty_breaks(n = 5)) +
labs(fill = "Frequency",
title = "Vehicle crime",
subtitle = "March 2017",
caption = "Contains OS data © Crown copyright and database right (2017)") +
theme_void()
library(plotly)
p <- ggplot(bdy_pts) +
geom_sf(aes(fill = n, text = paste0(ward, "\n", "Crimes: ", n))) +
scale_fill_gradientn('Frequency', colours=RColorBrewer::brewer.pal(5,"RdPu"),
breaks = scales::pretty_breaks(n = 5)) +
labs(fill = "Frequency",
title = "Vehicle crime",
subtitle = "March 2017",
caption = "Contains OS data © Crown copyright and database right (2017)") +
theme_void() +
coord_fixed(1.3)
ggplotly(p, tooltip = "text")
library(leaflet)
pal <- colorBin("RdPu", domain = bdy_pts$n, bins = 5, pretty = TRUE)
leaflet(bdy_pts) %>%
addTiles(urlTemplate = "http://{s}.tiles.wmflabs.org/bw-mapnik/{z}/{x}/{y}.png",
attribution = '© <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>, <a href="https://www.ons.gov.uk/methodology/geography/licences">Contains OS data © Crown copyright and database right (2017)</a>') %>%
addPolygons(fillColor = ~pal(n), fillOpacity = 0.8,
weight = 1, opacity = 1, color = "black",
label = ~as.character(ward)) %>%
addLegend(pal = pal, values = ~n, opacity = 0.7,
title = 'Vehicle crime (March 2017)', position = "bottomleft") %>%
addMiniMap(tiles = providers$CartoDB.Positron, toggleDisplay = TRUE)